
#include <cusp/blas.h>
#include <thrust/transform_reduce.h>

namespace cusp
{
namespace detail
{

template<typename ArrayType>
typename ArrayType::value_type frobenius_norm(ArrayType &A, cusp::array2d_format)
{
	return cusp::blas::nrm2(A.values);
}

template<typename MatrixType>
typename MatrixType::value_type frobenius_norm(MatrixType &A, cusp::hyb_format)
{
	typedef typename MatrixType::value_type ValueType;

	return sqrt(
		thrust::transform_reduce(A.ell.values.values.begin(), A.ell.values.values.end(),
			cusp::blas::detail::square<ValueType>(),
			ValueType(0),
			thrust::plus<ValueType>())
		+
		thrust::transform_reduce(A.coo.values.begin(), A.coo.values.end(),
			cusp::blas::detail::square<ValueType>(),
			ValueType(0),
			thrust::plus<ValueType>())
		);
}

template<typename MatrixType>
typename MatrixType::value_type frobenius_norm(MatrixType &A, cusp::coo_format)
{
	return cusp::blas::nrm2(A.values);
}

template<typename MatrixType>
typename MatrixType::value_type frobenius_norm(MatrixType &A, cusp::csr_format)
{
	return cusp::blas::nrm2(A.values);
}

template<typename MatrixType>
typename MatrixType::value_type frobenius_norm(MatrixType &A, cusp::sparse_format)
{
	return cusp::blas::nrm2(A.values.values);
}


} // end namespace detail


/////////////////
// Entry Point //
/////////////////

template<typename MatrixType>
typename MatrixType::value_type frobenius_norm(MatrixType &A)
{
	return detail::frobenius_norm(A, typename MatrixType::format());
}


} // end namespace cusp

